for j=1:length(ss)
s=ss(j);
betas=beta0+beta1*s; 
lambdas=max(lambda0*(1-zeta0*s^zeta1),0);

zetas=max(eta0*(1-s^eta1),0);
zetan(j)=zetas;
ix=ixx(j);
x=xx(j);
phii=ix-etaI*ix^2/2-deltaI;
 gx(j)=phii-lambdas/(betas+1)-zetan(j)*(qx(j)-qxB(j))/qx(j);
end

[xx(flagstar),ss(flagstar), ixx(flagstar), qx(flagstar),  cx(flagstar),gx(flagstar), rx(flagstar), rpx(flagstar), bx(flagstar)]